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Abstract 

An analysis of arrival directions of extensive air showers (EAS) registered with 
the EAS MSU and EAS-1000 prototype arrays has revealed a region of excessive 
flux of PeV cosmic rays in the direction toward pulsars PSR J1840+5640 and LAT 
PSR J1836+5925 at significance level up to 4.5o\ The first of the pulsars was 
discovered almost 30 years ago and is a well-studied old radio-pulsar located at 
the distance of 1.7 pc from the Solar system. The second pulsar belongs to a new 
type of pulsars, discovered by the space gamma-ray observatory Fermi, pulsations 
of which are not observed in optical and radio wavelengths but only in the gamma- 
ray range of energies (gamma-ray-only pulsars). In our opinion, the existence of 
the region of excessive flux of cosmic rays registered with two different arrays 
provides a strong evidence that isolated pulsars can give a noticeable contribution 
to the flux of Galactic cosmic rays in the PeV energy range. 

PACS codes: 98.70.Sa Cosmic rays (including sources, origin, acceleration, and in- 
teractions), 96.50.sd Extensive air showers, 97.60.Gb Pulsars, 98.70.Rz gamma-ray 
sources 



1 Introduction 

The problem of the nature of the knee in the cosmic ray (CR) energy spectrum at 
3 PeV and, in a more generic aspect, of the origin of CRs in the PeV energy range 
stays open for many years now. Attempts to collate arrival directions of CRs with 
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coordinates of their possible astrophysical sources is one of the routes of investigations 
on the subject. Another direction, which is being actively developed, is observation of 
gamma rays from such sources. 

In both approaches, the main subject of investigation are Galactic supernova rem- 
nants (SNRs) since the model of diffusive shock acceleration at the outer front of 
expanding SNRs has become widely spread, see, e.g., a review by Hillas 2005. This 
does not imply though that there are no other astrophysical sources in the Galaxy 
that are able to accelerate hadrons to PeV and higher energies. In particular, the 
first results demonstrating that pulsars can act as effective CR accelerators appeared 
soon after the discovery of pulsars and their identification with neutron stars, see, e.g., 
Gunn & Ostriker 1969. After this, interest to pulsars as possible sources of CRs with 
energies > 10 14 eV up to the highest energies did not vanish, see, e.g., Blasi et al. 2000; 
Giller & Lipski 2002; Bhadra 2003; Erlykin & Wolfendale 2004; Bhadra 2006. Models 
considered included acceleration both in pulsar wind nebulae (e.g., the Crab nebula) 
and by isolated pulsars. It was demonstrated in particular that the Geminga pulsar is 
a possible candidate for being the single source of the knee at 3 PeV (Bhadra 2003). 

In our previous works, we have already reported on the number of regions of exces- 
sive flux (REFs) found during an analysis of arrival directions of EAS registered with 
the EAS MSU array and the EAS-1000 prototype array ("PRO-1000") and generated 
by CRs with energies of around PeV, (Kulikov k, Zotov 2004; Zotov k, Kulikov 2004, 
2007, 2009). It was demonstrated that there are no Galactic SNRs in the vicinity 
of some of the REFs, though a considerable number of regions have isolated pulsars 
nearby. In Zotov Sz Kulikov 2009, we pointed out that there are two overlapping REFs 
in the vicinity of a point with equatorial coordinates a = 285°, 5 = 58° in the EAS 
MSU and the PRO-1000 data sets. Taking into account possible errors in determi- 
nation the boundaries of the REFs, one could consider them as one and the same 
region with the pulsar PSR J1840+5640 inside (coordinates of the pulsar a ~ 280.19°, 
5 ~ 56.68°). There are neither supernova remnants or SNR candidates inside or close 
to the region, nor open clusters or OB-associations, which are sometimes considered as 
possible sources of Galactic cosmic rays. PSR J1840+5640 was discovered back in 1980 
(Shitov et al. 1980) and is a well-studied old radio-pulsar of the age of 17.5 x 10 6 years 
located at 1.7 kpc from the Solar system, see, e.g., Arzoumanian et al. 1994; Hobbs 
et al. 2004. 

This region attracted our attention once again after there appeared a list of the 
brightest gamma-ray sources in the energy range from 100 MeV to 100 GeV registered 
during the first three month of operation of the Fermi Gamma-ray Space Telescope 
(Abdo et al. 2009a). It happened that there is a pulsar (coordinates a « 279.06°, 
5 59.46°) of a new type near the REF, a so called gamma-ray-only pulsar, i.e., a 
pulsar that only "blinks" in gamma rays. This pulsar as an intensive gamma-ray source 
was discovered by EGRET in 1991-1993 but the nature of the source remained unclear 
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(Nolan et al. 1994). It occurred to be the brightest unidentified gamma-ray source in 
the energy ranges > 100 MeV (significance level of 19. Oct) and > 1 GeV (14.2<r) located 
far from the Galactic plane (Galactic latitude b ps 25°), see Hartman et al. 1999; Lamb 
et al. 1997. 

Numerous attempts to elucidate the nature of the object named as 3EG J1835+5918 
and to find a definite astrophysical counterpart did not succeed, see Mirabal et al. 2000; 
Reimer et al. 2001; Halpern et al. 2002; Fegan et al. 2005; Halpern et al. 2007, and 
references therein. After intensive studies of its error box in optical, X-ray, and radio 
wave bands it was identified as an isolated neutron star correlating with a faint X-ray 
source RX J1836. 2+5925, discovered by ROSAT. It was suggested that the object is 
a Geminga-like pulsar, which is not observed in optical and radio wave bands, and an 
upper limit of the distance to the object was obtained (< 800 pc). All attempts to 
find pulsations in the X-ray and gamma-ray energy ranges failed. Pulsations of 3EG 
J1835+5918 were not observed during the first year of observations by the AGILE 
instrument either, though its gamma-ray emission was confirmed at 15. 6a significance 
level (Bulgarelli et al. 2008; Pittori et al. 2009). Only gamma-ray observations by the 
Fermi LAT lead to unequivocal identification of 3EG J1835+5918 as a pulsar, named 
LAT PSR J1836+5925. 

Since the search for gamma-ray sources had become a natural part of the search 
for PeV cosmic ray sources (Rowell et al. 2005), we performed a more detailed analysis 
of the REF located in the direction to both pulsars. 

2 Experimental Data 

Similar to our previous studies of EAS arrival directions, the main data sets consist 
of 513 602 showers registered with the EAS MSU array in 1984-1990 and 1 342 340 
showers of the PRO-1000 array registered in 1997-1999. Technical description of the 
arrays were given by Vernov et al. 1979, and Fomin et al. 1999, respectively. All EAS 
selected for the analysis satisfy a number of quality criteria and have zenith angles 
9 < 45.7°. Showers registered with the two arrays have different number of charged 
particles N c in a typical event. For the EAS MSU array, the median value of N e is of 
the order of 1.6 x 10 5 , while that for the PRO-1000 array equals 3.7 x 10 4 . Assuming 
that the majority of the EAS was generated by protons, one can estimate their energy 
as being at around 1 PeV. Modern models of hadronic interactions and data on the 
chemical composition of cosmic rays allow one to estimate energy of a primary particle 
that caused an EAS. The given values of N e correspond to the energies of a primary 
proton Eq ss 1.7 x 10 15 eV and Eq ps 4.6 x 10 14 eV respectively, with an error of about 
10%-20%. 

The EAS MSU array allowed one to determine arrival directions of EAS with a 



3 



better accuracy than the PRO-1000 array but the mean error is estimated to be of the 
order of 3° in both cases. 

3 Method of the Investigation 

The investigation is based on the method by Alexandreas et al. (Alexandreas et al. 
1991), which has been used multiple times for the analysis of arrival directions of 
CRs by different research groups. The idea of the method is as follows. To every 
shower in the experimental data set, arrival time of another shower is assigned in a 
pseudo-random way. After this, new equatorial coordinates (a, 5) are calculated for 
the "mixed" data set thus providing a "mixed" map of arrival directions. The mixed 
map differs from the original one but has the same distribution in declination 5. In 
order to compare both maps, one divides them into sufficiently small "basic" cells. A 
measure of difference between any two regions (cells) of the two maps located within 
the same boundaries is defined as 

S = (N re!il -N hg )/y/K~ g , 

where N Te3u \ and N^ g are the number of EAS inside a cell in the real and "mixed" 
(background) maps correspondingly. The mixing of the real map is performed multiple 
times, and the mixed maps are averaged then in order to reduce the dependence of 
the result on the choice of arrival times. The method is based on an assumption 
that the resulting mean "background" map has most of the properties of an isotropic 
background, and presents the distribution of arrival directions of cosmic rays that 
would be registered with the array in case there is no anisotropy. Thus, deviations 
of the real map from the background one may be assigned to a kind of anisotropy 
of arrival directions of EAS registered at the array. As a rule, selection of regions of 
excessive flux is performed basing on the condition S > 3. 

Similar to the paper Zotov & Kulikov 2009, basic cells are of the size 0.5° x 0.5°. 
The number of cycles of mixing was increased from 1000 to 10,000 in order to further 
improve the quality of the background map. Due to this, the difference between the 
basic cells of any two consecutive averaged maps obtained close to the end of the mixing 
reduced approximately by an order of magnitude and is ~ 1-2 x 10~ 3 for both data sets. 

Cells of excessive flux (CEFsfl of CRs were searched for in the following way. 
Adjacent basic strips in 5 (each 0.5° wide) were joined into strips of width A5 = 
3° . . . 30° with step equal to 0.5°. Each wide strip was then divided into adjacent cells 
of some fixed width Aa. After this, we calculated the number of EAS inside each of 

1 We use the word "cell" for regions selected during the initial analysis and preserve the term "region 
of excessive flux" for regions formed of multiple cells. Thus, contrary to cells, REFs are not necessarily 
quadrangles. 
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these cells for both experimental (N rea \) and background (iV bg ) maps. For each pair of 
cells, S was calculated. A cell was considered as a possible CEF if S > 3. Every wide 
strip was divided into cells with a shift of 0.5° in a until all possible locations on the 
grid were covered. 

Contrary to Zotov & Kulikov 2009, we did not restrict the search to cells such 
that Act = AS/ cos 5, where 5 is the mean value of 5 for the current strip, and Aa is 
rounded to the nearest half-integer number, but allowed Aa to deviate from that value 
by Aa/6. As a result, additional cells were considered. For example, besides cells of 
the size Aa x A5 = 3° x 3°, we also considered cells of the sizes 2.5° x 3° and 3.5° x 3°. 
Finally, only cells with more than 10 EAS inside were analyzed. 

The method of Alexandreas et al. does not provide a direct answer to the question 
about the chance probability of appearance of a CEF. It is implicitly assumed that the 
deviation of AT real from Abg has a Gaussian distribution, and thus one is expected to 
calculate the corresponding probability basing on the value of S, which is an estimate 
of the standard deviation of a sample and thus acts as a significance level. Hence, it is 
assumed the chance probability of a CEF to appear is less than 1 — 0.9973 providing 
it was selected at significance level S > 3. 

In order to estimate the chance probability of appearance of a CEF basing on 
the number of EAS inside, we introduced the following simple method based on the 
binomial distribution (Zotov & Kulikov 2009). Let a shower axis getting inside the CEF 
be a success. The number of trials equals the number of showers N in the data set 
under consideration, and an estimate of success (for a fixed region) equals p = N^g/N, 
where N^ g is the number of showers in the cell of the background map. The assumption 
is based on the fact that in the method of Alexandreas et al. , Ab g is considered as an 
expected number of showers in a cell. Obviously, the chance probability of finding 
exactly A rca i EAS in a cell (or region) equals 

Piy = A rcal ) = C(N, N vcal )p N ^(l - p) N - N -\ 

where v is a random variable equal to the number of successes in the binomial model, 
and C(A, A rca i) is the corresponding binomial coefficient. Obviously, it is more in- 
teresting to consider a probability that there are at most iV rea i showers in the cell 
P/v = P{y < A rca i), or its adjacent value 1 — Pjv- The analysis performed and the data 
presented below demonstrate that values of Pn correlate well with the values of chance 
probabilities calculated on the basis of significance level S. 

The fact that S > 3 for a given cell does not imply it is valid for its sub-cells. 
More than this, since the method of Alexandreas et al. does not anyhow restrict sizes 
and shapes of regions to be considered, the non-uniformity of the distribution of data 
w.r.t. declination 5 can result in a selection of a CEF such that a great value of S is 
obtained due to an excess of EAS in just one of its sub-cells. A situation such that 
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S < 3 for a number of adjacent cells but S > 3 for a cell (region) combined of them 
is possible. In order to improve the robustness of selection of CEFs, we have tried a 
number of additional quantities. One of the most efficient of them is the probability 
Pbc = P(£, < ^bc)i calculated from the following binomial model. Let £ be a random 
variable equal to the number of basic cells of the given cell with an excess of EAS 
over the background values, and N& is the corresponding experimental value. The 
number of trials equals the number of basic cells in the CEF. It is natural to assume 
the probability of success to be equal to 1/2. In the results presented below, all CEFs 
satisfy a condition P\^ c > 0.9545, which corresponds to the significance level of 2a for 
the Gaussian distribution. We thus reduced the chance that a CEF is selected solely 
due to the non-uniformity of the EAS distribution w.r.t. 5. 

4 Results 

The region of excessive flux of CRs in the direction toward pulsars PSR J1840+5640 
and LAT PSR J1836+5925 found by the analysis of the main EAS MSU data set is 
shown in the upper panel of Figure 1. The REF consists of 99 CEFs with S > 3.0 
located within the boundaries a = 277.5° . . . 289°, S = 55.5° . . . 61°. Some parameters 
of the CEFs the region consists of are presented in Table 1. Pbc = 0.989308 for a CEF 
with the highest value of S. For this CEF, an excess of the experimental flux over 
the background one takes place in 53 of 84 basic cells that form the CEF. The set 
of showers within the REF has a slightly greater median value of N c than the whole 
data set. For example, median N c « 1.7 x 10 5 for 1508 EAS with a = 277.5° . . . 289°, 

5 = 56° ... 60°, and median N c » 1.8 x 10 5 for 856 EAS within the region covered by 
CEFs with S > 4.3. 

The result obtained for the whole EAS MSU data set is confirmed in the case 
one selects subsets with an improved accuracy of determination of arrival directions. 
The most highly pronounced CEF for a set of showers with N e > 10 5 (343 141 EAS) 
has S = 3.951, P N = 0.999947, P hc = 0.993082 and contains 620 EAS within the 
boundaries a = 280.5° . . . 287°, 5 = 54.5° . . . 58.5° with the background level equal to 
529.1 EAS. For this CEF, an excess of the experimental flux over the background one 
takes place in 65 of 104 basic cells that form the CEF. 

A similar result is obtained if one excludes showers that have arrival directions close 
to the vertical. In the set of showers with the zenith angle 6 > 6° (492 853 EAS), we 
found a CEF with S = 3.985, P N = 0.999952, P bc = 0.961815 within the boundaries 
a = 280° . . . 286°, 6 = 56° ... 59°. For the CEF, an excess of the experimental flux 
over the background one takes place in 44 of 72 basic cells. Thus, in all these cases 
we observe regions of excessive flux with the probability to appear by chance that can 
be estimated as 7.8 x 10~ 5 -7.5 x 10" 6 if one uses the maximum values of S, or as 
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Figure 1: Upper panel: the region of excessive flux of CRs found in the main data set 
of the EAS MSU array. Lower panel: the region of excessive flux of CRs found in the 
PRO-1000 data set for showers with N c > 3 x 10 4 . 



5.3 x 10 -5 -6 x 10~ 6 in case one uses the corresponding values of Pjv- 

A region of excessive flux of CRs in the direction toward pulsars PSR J1840+5640 
and LAT PSR J 1836+5925 was also found in the full data set of the PRO-1000 array 
(with the significance level up to S = 3.573 and Pn = 0.999800) and in its subsets 
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Table 1: Parameters of cells of excessive flux that form the REF found in the EAS 
MSU data set (see Figure 1). Notation: S is the significance level for selecting CEFs, 
Nqef is the number of CEFs selected for a given value of S, a and S are the ranges 
of the corresponding equatorial coordinates, -/Veas is the range of the number of EAS 
registered in selected CEFs, Ab g is the corresponding range background events, min P/v 
is the minimum value of P/v for the CEFs selected. 
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similar to those described above. Namely, we found a CEF with S = 3.785, Pjv = 
0.999908 in a subset with 9 > 6° (1 267 089 EAS), and a CEF with S = 3.603 for a 
subset with N c > 3 x 10 4 (471 554 EAS). A REF found in the latter case is shown in the 
lower panel of Figure 1. Some parameters of CEFs that form the REF are presented 
in Table 2. The median value of N c for the CEF selected at the highest significance 
level equals 6.9 x 10 4 . This is slightly less than the corresponding value for the whole 
set of showers with N c > 3 x 10 4 equal to 7.2 x 10 4 . 



Table 2: Parameters of CEFs that form a region of excessive flux of CRs in the subset 
of the PRO-1000 data with N c > 3 x 10 4 . See notation in Table 2. 
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0.999811 



Thus, we have found cells of excessive flux in the experimental data of the PRO- 
1000 array with the probabilities of appearing by chance that can be estimated as 
1.5-3.5 xl0~ 4 if we use the maximum values of the significance level S, or as l-2xl0~ 4 
in case estimates are based on the corresponding values of Pjv- 

There is a possibility that coordinates of the two pulsars coincide with the position 
of the REF just by chance. We think a probability of such a coincidence is small 
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because the majority of known pulsars are located in the vicinity of the Galactic plane 
while the pulsars under consideration have galactic latitude b ~ 25°. In particular, less 
than 20% of more than 330 pulsars known in the region S > 10° have b > 20°, and only 
26 of them are located at distances less than 2 kpc from the Solar system (Manchester 
et al., 2005). 

Flux of primary cosmic rays provides additional information about the REF found. 
One must take into account though that an area covered by an EAS array does not 
coincide with the so called "effective" area. The EAS MSU array covered an area of 
approximately 0.5 km 2 but a circle of effective registration was defined for the proba- 
bility of registering a shower > 0.95. It has a radius equal to 20 m for showers with 
N c > 2.5x 10 5 (Eq > 2.5 PeV), 30 m for showers with N e > 4x 10 5 (E > 3.7 PeV), etc. 
An area of effective registration for the PRO-1000 array was defined less accurately 
thus we only provide data for the EAS MSU array. The flux of primary cosmic rays in 
the REF equals ps 4.4 x 10" 13 cm" 2 s" 1 for showers with A" e =2.5-4xl0 5 (237 EAS), 
and « 2.3 x 10~ 13 cm" 2 s" 1 for showers with A" e =4xl0 5 -10 6 (273 EAS.) The number 
of showers in the REF with greater values of N e is less than 100, thus we do not take 
it as statistically sufficient. 

5 Discussion 

In spite of the fact that the boundaries of the regions of excessive flux found in the 
data sets obtained with two different experimental arrays are shifted with respect to 
each other, we think they point to one and the same region. Making this conclusion, 
we take into account the discreetness of the grid and possible errors in determination 
of arrival directions of EAS. 

As we have already mentioned in the introduction, there are no other Galactic 
objects usually considered as possible sources of PeV cosmic rays in the vicinity of the 
REF described above. The list of objects checked includes supernova remnants (Green 
2009), including pulsar wind nebulae (Roberts 2005), open clusters (Dias et al. 2002), 
OB-associations, and SNR candidates (SIMBAD database). Hence, a natural question 
is whether pulsars PSR J1840+5640 and LAT PSR J1836+5925 are able to accelerate 
protons and/or heavier nuclei up to energies of the order of PeV. 

To estimate the maximum energy of a particle accelerated in the wind near the 
light cylinder of a pulsar, one can use the following expression (Blasi et al. 2000): 
£ max = 0.34 Z BQ, 2 eV, where Z is the charge number of the particle, B is the strength 
of the surface magnetic field, G, Q is the angular velocity of a pulsar, rad s _1 , and the 
radius of the pulsar equals 10 6 cm. The typical energy of the accelerated particle equals 
approximately a half of E max . 

For PSR J1840+5640, B » 1.59 x 10 12 G, fi as 3.80 rad s" 1 (barycentric period 
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P ps 1.653 s, see Manchester et al. 2005). Thus, the maximum energy of a proton 
accelerated by this pulsar equals approximately 7.8 x 10 12 eV. An iron nuclei can be 
accelerated up to 2 x 10 14 eV. It follows that PSR J1840+5640 can contribute to the 
formation of the REF only in case there is a considerable fraction of heavy nuclei in 
the flux of primary CRs from the corresponding direction. 

For LAT PSR J1836+5925, B « 0.5 x 10 12 G, n » 36.26 rad s^ 1 (Abdo et al. 
2009b). Thus, the maximum energy of a proton accelerated by this pulsar is of the 
order of 2.2 x 10 14 eV. For an iron nuclei, E mSbX « 5.8 x 10 15 eV. Hence, LAT PSR 
J1836+5925 is a possible candidate for being an accelerator of charged particles in the 
PeV energy range. 

Let us assume basing on the estimates of E nmx that the REF described above was 
formed due to the contribution of one or both pulsars. Then there arises a question of 
how was the arrival direction of the corresponding EAS kept so close to the direction 
to these pulsars. It is difficult to give a definite answer. On the one hand, we cannot 
totally exclude a possibility that a part of showers that formed the REF was initiated 
by photons but not hadrons. The PRO-1000 array did not allow one to distinguish 
between gamma-ray and hadronic showers but the EAS MSU array provided such an 
opportunity because it had a muon detector. Since photonic showers contain very few 
muons in comparison with hadronic ones, their fraction in the whole data set could 
be estimated. It was shown by Khristiansen et al. (1975) that the fraction of photon- 
initiated showers in the data set obtained with the EAS MSU array in 1960s was of the 
order of 10 -3 . A similar estimate was obtained basing on the results of the CASA-MIA 
experiment, see Chantell et al. 1997. Therefore, we have no ground to assume that the 
observed excess of CRs is formed due to photons. 

Another possibility is a negligible influence of the Galactic magnetic field (GMF) 
on particle trajectories, or even magnetic lensing. It is difficult to say anything definite 
about the transverse component of the GMF along the line of sight to the pulsars but 
provide the usual estimates of the mean values of the field in the vicinity of the Solar 
system. It is remarkable though that there exists a local magnetic field going toward 
the Galactic longitude I 90°, close to that of both pulsars, with a strength about 2- 
3 fiG (Manchester 1974). Known values of the rotation measure RM and the dispersion 
measure DM of the pulsar PSR J1840+5640 allow one to estimate the average value 
of the GMF along the line of site to the pulsar (see, e.g., Han et al. 2006): (Bii) = 
1.232 RM/DM fiG. For PSR J1840+5640, RM = -3 rad m~ 2 , DM = 26.698 cm" 3 pc 
(Manchester et al. 2005). Thus, (B\\) rj -0.138 fiG. 

One might ask, why has not been the REF found in experimental data of other 
arrays that studied PeV cosmic rays. Actually, one can find a number of reports about 
small-scale anisotropy (regions of excessive flux of CRs) in this energy range, see, e.g., 
Benko et al. 2004, and Sun Sz Sun 1997. Some of the regions presented in these works 
have coordinates close to the other regions found by us, but mostly the results do not 
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agree. 

On the other hand, a search for point-like sources of PeV CRs with the data of 
CASA-MIA and KASCADE arrays gave a negative result, see McKay et al. 1993, 
Antoni et al. 2004. A method close to that of Alexandreas et al. was used for the 
analysis of the CASA-MIA data. A comparison of the experimental and background 
fluxes was performed for regions of the size 0.8° x 0.8°/ cos 5. An existence of at least 
four adjacent regions with an excess of the number of registered showers over the 
background at > 2.8a was demanded. The KASCADE data set was studied with the 
original method by Alexandreas et al. Only regions of the size 0.5° x 0.5°/ cos 5 were 
analysed. It is possible that in both cases the choice of such small regions was the main 
reason of the negative result: the fraction of EAS generated by PeV photons is not 
sufficient for providing a considerable excess of the observed flux over the background, 
and it is unlikely that the arrival direction of a charged PeV particle is preserved 
with such a high precision. In our opinion, it would be very interesting to revisit the 
data of the CASA-MIA and KASCADE experiments and to study regions of the size 
> 3° x 3°/ cos 5. 

To conclude, we cannot unequivocally claim that the region of excessive flux of PeV 
cosmic rays described above was formed due to the contribution of any or both of the 
pulsars PSR J1840+5640 and LAT PSR J1836+5925. Still, we think the existence of 
the REF in the direction toward these pulsars is worth paying attention as a possible 
sign that isolated pulsars provide a more noticeable contribution to the flux of Galactic 
cosmic rays than it is sometimes assumed. This point of view is supported by the fact 
that there are four more REFs in the EAS MSU data set selected with significance 
level S > 4 that can be associated with pulsars which do not have known SNRs nearby. 
These REFs will be the subject of a separate article. 
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